Part 1: compute relative tidal range with respect to waves
Contents
import pathlib
import sys
# make coastpy library importable by appending root to path
cwd = pathlib.Path().resolve()
proj_dir = cwd.parent.parent
sys.path.append(str(proj_dir / "src"))
# definitly not a best practice, but a workaround to avoid many warnings caused by shapely 2.0 release
import io
import socket
import warnings
import cartopy.crs as ccrs
import geopandas as gpd
import geoviews as gv
import geoviews.tile_sources as gts
import geoviews.tile_sources as gvts
import holoviews as hv
import hvplot.pandas
import hvplot.xarray
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import panel as pn
import requests
import rioxarray
import shapely
import xarray as xr
from tqdm import tqdm
from coastpy.geometries import geo_bbox
warnings.filterwarnings("ignore", category=RuntimeWarning)
warnings.filterwarnings("ignore", category=UserWarning)
/var/folders/p8/qxnqn2ns5kbczrp91kx4_1nm0000gn/T/ipykernel_14066/1733808755.py:7: UserWarning: Shapely 2.0 is installed, but because PyGEOS is also installed, GeoPandas will still use PyGEOS by default for now. To force to use and test Shapely 2.0, you have to set the environment variable USE_PYGEOS=0. You can do this before starting the Python process, or in your code before importing geopandas:
import os
os.environ['USE_PYGEOS'] = '0'
import geopandas
In a future release, GeoPandas will switch to using Shapely by default. If you are using PyGEOS directly (calling PyGEOS functions on geometries from GeoPandas), this will then stop working and you are encouraged to migrate from PyGEOS to Shapely 2.0 (https://shapely.readthedocs.io/en/latest/migration_pygeos.html).
import geopandas as gpd
DATA_DIR = proj_dir / "data"
coastsat_metadata_pacific_fp = DATA_DIR / "03_CoastSat_metadata_layer_Pacific.geojson"
metadata = gpd.read_file(coastsat_metadata_pacific_fp)
for c in ["Tier1Images", "Tier2Images", "S2Images", "pk"]:
metadata[c] = metadata[c].astype("i8")
for c in ["TideRange", "WaveHeight", "RelativeTideRange"]:
metadata[c] = metadata[c].astype("f8")
metadata.head()
| Tier1Images | Tier2Images | S2Images | TideRange | WaveHeight | RelativeTideRange | pk | geometry | |
|---|---|---|---|---|---|---|---|---|
| 0 | 410 | 113 | 207 | 1.65 | 1.22 | 1.35 | 0 | POINT (180.00000 -16.15293) |
| 1 | 256 | 33 | 214 | 1.79 | 0.60 | 2.97 | 1 | POINT (179.16992 -16.41845) |
| 2 | 255 | 56 | 415 | 1.86 | 0.81 | 2.29 | 2 | POINT (178.51402 -16.78987) |
| 3 | 301 | 78 | 214 | 1.47 | 0.60 | 2.43 | 3 | POINT (179.18860 -16.72214) |
| 4 | 495 | 176 | 214 | 1.51 | 1.01 | 1.49 | 4 | POINT (179.89177 -16.65793) |
Part 1: compute relative tidal range with respect to waves#
todo: check with Kilian which waves heights he used mean_tidal_range / mean_wave_height
metadata["RelativeTideRange"] = metadata["TideRange"] / metadata["WaveHeight"]
# metadata["RelativeTideRange"] = ... # your code to compute tidal range over the wave height
pn.extension()
title_bar = pn.Row(
pn.pane.Markdown(
"## Tide versus wave dominated coasts",
style={"color": "black"},
width=500,
sizing_mode="fixed",
margin=(10, 5, 10, 15),
),
pn.Spacer(),
)
w_threshold = pn.widgets.FloatSlider(
name="Threshold", start=0.1, end=10.0, step=0.1, value=1
)
@pn.depends(w_threshold.param.value)
def plot(threshold):
wave_dominated = metadata[metadata["RelativeTideRange"] < threshold].copy()
tide_dominated = metadata[metadata["RelativeTideRange"] > threshold].copy()
return (
gts.EsriImagery()
* wave_dominated.hvplot(geo=True, color="blue", label="Wave dominated")
* tide_dominated.hvplot(geo=True, color="green", label="Tide dominated")
).opts(legend_position="bottom_right")
app = pn.Column(
title_bar,
pn.Row(w_threshold),
pn.Row(plot),
)
app
from ipyleaflet import Map, basemaps
m = Map(basemap=basemaps.Esri.WorldImagery, scroll_wheel_zoom=True)
m.center = 37.768137, -122.511066
m.zoom = 16
m.layout.height = "800px"
m
bbox = geo_bbox(m.west, m.south, m.east, m.north)
bbox.explore()
Make this Notebook Trusted to load map: File -> Trust Notebook
transects = gpd.read_file(DATA_DIR / "03_CoastSat_transect_layer.geojson")
transects_roi = gpd.sjoin(transects, bbox)
transects_roi.explore()
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[10], line 1
----> 1 transects_roi.explore()
File ~/mambaforge/envs/codebook/lib/python3.10/site-packages/geopandas/geodataframe.py:1981, in GeoDataFrame.explore(self, *args, **kwargs)
1978 @doc(_explore)
1979 def explore(self, *args, **kwargs):
1980 """Interactive map based on folium/leaflet.js"""
-> 1981 return _explore(self, *args, **kwargs)
File ~/mambaforge/envs/codebook/lib/python3.10/site-packages/geopandas/explore.py:366, in _explore(df, column, cmap, color, m, tiles, attr, tooltip, popup, highlight, categorical, legend, scheme, k, vmin, vmax, width, height, categories, classification_kwds, control_scale, marker_type, marker_kwds, style_kwds, highlight_kwds, missing_kwds, tooltip_kwds, popup_kwds, legend_kwds, map_kwds, **kwargs)
363 map_kwds["max_zoom"] = tiles.get("max_zoom", 18)
364 tiles = tiles.build_url(scale_factor="{r}")
--> 366 m = folium.Map(
367 location=location,
368 control_scale=control_scale,
369 tiles=tiles,
370 attr=attr,
371 width=width,
372 height=height,
373 **map_kwds,
374 )
376 # fit bounds to get a proper zoom level
377 if fit:
File ~/mambaforge/envs/codebook/lib/python3.10/site-packages/folium/folium.py:274, in Map.__init__(self, location, width, height, left, top, position, tiles, attr, min_zoom, max_zoom, zoom_start, min_lat, max_lat, min_lon, max_lon, max_bounds, crs, control_scale, prefer_canvas, no_touch, disable_3d, png_enabled, zoom_control, **kwargs)
272 zoom_start = 1
273 else:
--> 274 self.location = validate_location(location)
276 Figure().add_child(self)
278 # Map Size Parameters.
File ~/mambaforge/envs/codebook/lib/python3.10/site-packages/folium/utilities.py:78, in validate_location(location)
71 raise ValueError(
72 "Location should consist of two numerical values, "
73 "but {!r} of type {} is not convertible to float.".format(
74 coord, type(coord)
75 )
76 )
77 if math.isnan(float(coord)):
---> 78 raise ValueError("Location values cannot contain NaNs.")
79 return [float(x) for x in coords]
ValueError: Location values cannot contain NaNs.
TIMEOUT_SEC = 5 # default timeount in seconds
# Get the transect id's by saving the transects in our region of interest as a list.
transect_ids = transects_roi.TransectId.to_list()
# Here we download the data per transect, but if it takes too long we default to Ocean beach, California.
try:
data = []
for transect_id in tqdm(transect_ids):
url = url = f"http://coastsat.wrl.unsw.edu.au/time-series/{transect_id}/"
# Pandas.read_csv() can handle url's but it doesn't accept a timeout argument,
# so we download the data via requests.
r = requests.get(url, timeout=TIMEOUT_SEC).content
s = pd.read_csv(
io.StringIO(r.decode("utf-8")),
header=None,
names=["date", "shoreline_position"],
)
s["date"] = pd.to_datetime(s["date"])
s = s.set_index("date").replace({"None": np.nan}).astype("f8")
s.name = transect_id
data.append(s)
column_names = [s.name for s in data]
shorelines = pd.concat(data, axis=1)
shorelines.columns = column_names
except:
# If it takes too long to get the data from the server, default to Ocean Beach, California.
# First set the map to that area, so that we have the correct bbox.
m.center = 37.768137, -122.511066
m.zoom = 16
m.layout.height = "800px"
bbox = geo_bbox(m.west, m.south, m.east, m.north)
transects_roi = gpd.sjoin(transects, bbox)
shorelines_fp = DATA_DIR / "03_shorelines_ocean_beach.parquet"
shorelines = pd.read_parquet(shorelines_fp)
0%| | 0/8 [00:00<?, ?it/s]
shorelines.head()
| usa_CA_0215-0003 | usa_CA_0215-0004 | usa_CA_0215-0005 | usa_CA_0215-0006 | usa_CA_0215-0007 | usa_CA_0215-0008 | usa_CA_0215-0009 | usa_CA_0215-0010 | usa_CA_0215-0011 | usa_CA_0215-0012 | usa_CA_0215-0013 | usa_CA_0215-0014 | usa_CA_0215-0015 | usa_CA_0215-0016 | usa_CA_0215-0017 | usa_CA_0215-0018 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| date | ||||||||||||||||
| 1984-05-02 18:13:15 | 71.4 | 87.6 | 103.9 | 109.0 | 115.4 | NaN | 120.1 | 124.8 | NaN | NaN | NaN | NaN | 148.7 | 146.1 | 131.7 | NaN |
| 1984-06-03 18:14:10 | 74.4 | 70.1 | 98.5 | 116.8 | 117.9 | 141.0 | 127.4 | 134.9 | NaN | 166.6 | NaN | NaN | NaN | 139.3 | 133.6 | NaN |
| 1984-06-10 18:20:29 | 62.9 | 71.4 | 90.4 | 118.9 | 111.7 | 137.4 | 140.6 | 126.9 | 137.5 | 158.1 | 150.3 | 150.1 | 155.8 | 150.8 | 142.7 | NaN |
| 1984-06-26 18:20:42 | 60.7 | 74.3 | 91.6 | 112.9 | 127.7 | 133.8 | 152.4 | 156.0 | 157.7 | 171.2 | 164.9 | 174.4 | 165.9 | 148.2 | 144.7 | NaN |
| 1984-09-07 18:16:09 | 77.3 | 91.1 | 110.3 | 128.9 | 141.6 | 152.9 | 166.0 | 173.3 | 173.5 | 173.8 | 174.8 | 184.2 | 195.9 | 202.3 | 190.7 | 185.3 |
check this for inspiration on how to clor the plot> https://foundations.projectpythia.org/core/xarray/enso-xarray.html#
pn.extension()
title_bar = pn.Row(
pn.pane.Markdown(
"## Historical shoreline series",
style={"color": "black"},
width=500,
sizing_mode="fixed",
margin=(10, 5, 10, 15),
),
pn.Spacer(),
)
options = transects_roi.TransectId.to_list()
transect_slider = pn.widgets.Select(
name="Transect", options=options, value=np.random.choice(options)
)
@pn.depends(transect_slider.param.value)
def plot_transects(transect_id):
transect = transects_roi.loc[transects_roi["TransectId"] == transect_id].copy()
return (
gts.EsriImagery()
* transects_roi.hvplot(geo=True, color="blue")
* transect.hvplot(geo=True, color="red")
)
@pn.depends(transect_slider.param.value)
def plot_series(transect_id):
shoreline_selected = shorelines.hvplot(x="date", y=[transect_id], color="red")
shorelines_ = shorelines[
shorelines.columns.difference([transect_slider.value])
].hvplot(x="date", alpha=0.4)
return shorelines_ * shoreline_selected
app = pn.Column(
title_bar,
pn.Row(transect_slider),
pn.Row(pn.Column(plot_transects), plot_series),
)
app
mei_fp = DATA_DIR / "03_meiv2.data"
mei = (
pd.read_csv(mei_fp, delim_whitespace=True, index_col=0, header=None)
.reset_index()
.rename({0: "year"}, axis="columns")
.melt(id_vars="year", var_name="month")
)
mei["date"] = pd.to_datetime(mei["year"].astype(str) + "-" + mei["month"].astype(str))
mei = (
mei.drop(["year", "month"], axis="columns")
.set_index("date")
.sort_index()
.replace({-999: np.nan})
)
mei.hvplot(x="date")
mei_interp = mei.reindex(shorelines.index, method="nearest")
shorelines_ = shorelines.copy()
shorelines_.index = mei_interp.value.to_list()
shorelines_.hvplot(kind="scatter")
shorelines_ = shorelines.diff().copy()
shorelines_.index = mei_interp.value.to_list()
shorelines_.hvplot(kind="scatter", alpha=0.2, color="blue")